
#include "CUDA_interpolation.h"
#include <stdio.h>
#include "utility.h"

__device__ int linear_index(int x, int y, int dim){ return y * dim + x; }

__global__ void interpolation_kernel(float * d_coarse, float * d_fine, int fine_dim, int coarse_dim)
{
	//1-indexed coarse indices
	int coarse_x = blockIdx.x * blockDim.x + threadIdx.x + 1;
	int coarse_y = blockIdx.y * blockDim.y + threadIdx.y + 1;
	if(coarse_x <= coarse_dim && coarse_y <= coarse_dim)
	{
		float value = d_coarse[linear_index(coarse_x-1, coarse_y-1, coarse_dim)];
		float half_val = 0.5 * value;
		float quarter_val = 0.25 * value;
		//(2i, 2j)
		d_fine[linear_index(2*coarse_x-1, 2*coarse_y-1, fine_dim)] = value;

		//(2i-1, 2j)
		d_fine[linear_index( 2*coarse_x-2, 2*coarse_y-1, fine_dim )] += half_val;

		//(2i+1, 2j)
		d_fine[linear_index( 2*coarse_x, 2*coarse_y-1, fine_dim )] += half_val;

		//(2i, 2j-1)
		d_fine[linear_index(2*coarse_x-1, 2*coarse_y-2, fine_dim)] += half_val;

		//(2i, 2j+1)
		d_fine[linear_index(2*coarse_x-1, 2*coarse_y, fine_dim)] += half_val;

		//(2i-1, 2j-1)
		d_fine[linear_index(2*coarse_x-2, 2*coarse_y-2, fine_dim)] += quarter_val;

		//(2i-1, 2j+1)
		d_fine[linear_index(2*coarse_x-2, 2*coarse_y, fine_dim)] += quarter_val;

		//(2i+1, 2j-1)
		d_fine[linear_index(2*coarse_x, 2*coarse_y-2, fine_dim)] += quarter_val;

		//(2i+1, 2j+1)
		d_fine[linear_index(2*coarse_x, 2*coarse_y, fine_dim)] += quarter_val;

	}
}

void CUDA_Interpolation::operator()(CUDA_vec<float> & a, CUDA_vec<float> & b)
{
	int coarse_dim = (int)sqrt(a.getSize());
	int fine_dim = 2 * (coarse_dim+1) - 1;
	int BLOCK_SIZE = 16;
	int num_blocks_dim = (int)ceil(float(coarse_dim) / BLOCK_SIZE);
	b.set_constant(0);
//	printf("N_c: %d N: %d Block_size: %d\n",coarse_dim, fine_dim, num_blocks_dim);
	dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
	dim3 dimGrid(num_blocks_dim, num_blocks_dim);

//	cudaPrintfInit();
	interpolation_kernel<<<dimGrid, dimBlock>>>(a.raw_pointer(),b.raw_pointer(), fine_dim, coarse_dim);
//	cudaPrintfDisplay(stdout, true);
//	cudaPrintfEnd();
	//print_CUDA_err_msg();
	cudaThreadSynchronize();
}

